Numerical study of freezing in polydisperse colloidal suspensions

We have established the solid–fluid coexistence region for a system of polydisperse hard spheres with near-Gaussian diameter distributions, as a function of polydispersity. Significantly, we observe a ‘terminal’ polydispersity above which there can be no fluid–solid coexistence. At the terminus the polydispersity is 5.7% for the solid and 11.8% for the fluid while the volume fractions are 0.588 and 0.547 for the solid and fluid respectively. Substantial fractionation observed at high values of the polydispersity (>5%) implies that the ‘constrained eutectic’ assumption made in previous theoretical studies is not generally valid. Our results for the terminal polydispersity are consistent with experiments performed on polydisperse colloidal suspensions. Polydispersity in size is an unavoidable feature of colloidal systems, and thus it is of interest to understand the extent to which polydispersity can influence the behaviour of these systems. In this work we focus in particular on the crystallization properties of polydisperse suspensions. The hard-sphere model gives a very good description of colloidal dispersions that consist of non-charged spherical particles interacting via a steep steric repulsion [1]. In particular, such colloidal systems, if sufficiently monodisperse in size, are known to crystallize at densities very close to that predicted by a hard-sphere model [2]. However, it has been found that size polydispersity can greatly influence the location, and even the existence, of this transition. Given its simplicity and utility in describing nearly monodisperse suspensions, it is reasonable to extend the hard-sphere model in an examination of the effect of size polydispersity on the crystallization transition. The influence of polydispersity on the solid–fluid transition has been examined with theory [3–5], molecular simulation [6–8], and experiment [2]. A general conclusion from all of these studies is the existence of an upper limit on the size polydispersity, above which no crystallization can occur. We have chosen to call this the ‘terminal polydispersity’ [9] If the polydispersity is defined as the standard deviation of the size distribution divided by the mean, then the various studies indicate that the terminal polydispersity is in the range of 5 to 15 per cent. All of the studies based in theory or molecular simulation invoke some assumption regarding the distribution of particle diameters between the two phases. Some studies, but not all [4], assume that the diameter distribution is the same in both phases. In principle this is incorrect, as phase equilibria in mixtures generally involve some fractionation, with components of the mixture distributing themselves unevenly between the phases. The ‘constrained eutectic’ [5] approximation can be justified in the present 0953-8984/96/479627+05$19.50 c � 1996 IOP Publishing Ltd 9627 9628 P G Bolhuis and D A Kofke context if one focuses on nearly monodisperse systems, as much of the previous work has done. In doing so, however, one is left only to examine the influence of polydispersity on the location of the freezing point. As the presence and degree of fractionation is of some practical interest, it is worthwhile to attempt a more rigorous treatment that incorporates the effect. We have also found that the presence of size fractionation can explain some of the disparity in the values of the terminal polydispersity described above. This paper summarizes a recent study on the effect of size polydispersity on the freezing transition of the hard-sphere model [9]. We trace out the phase diagram for increasingly polydisperse systems, beginning with the well-studied [10, 11] limit of monodisperse hard spheres. Our approach is based on Monte Carlo simulation in a semi-grand ensemble, with the Gibbs–Duhem integration technique applied to trace out the coexistence lines [12]. Previous work [13, 14] has established the utility of a semi-grand ensemble in characterizing and simulating polydisperse mixtures. A simulation in the semi-grand ensemble has the total number of particles fixed, but the species identity of each particle is allowed to change, giving rise to a truly continuous distribution. Although the chemical potentials are imposed in a way similar to the grand canonical ensemble, insertion of particles is avoided, so the method is suitable for high densities and crystalline phases. For a system of N hard spheres with diameters σ distributed according to p(σ ), the differential form of the isobaric semi-grand canonical free energy Y is d(βY ) = H dβ + βV dP −N � p(σ )βδ[μ(σ )− μ(σ0)] dσ + βμ(σ0) dN. (1) Here, μ(σ ) is the chemical potential as a function of σ , σ0 is the diameter of an arbitrarily chosen reference component, and δ represents a functional differential. Also, H is the enthalpy, β = 1/kBT is the reciprocal temperature, P is the pressure and V the volume of the system. The (isobaric) semi-grand canonical potential Y is a function of the independent variables T , P and N and it is a functional of the chemical potential difference function �μ(σ ) ≡ μ(σ ) − μ(σ0). In a simulation these independent variables must be fixed while the thermodynamic conjugates H, V, μ(σ0) and p(σ ) are allowed to fluctuate. In nearly monodisperse systems, it can be expected that the form of the composition distribution p(σ ) will be much like that of the imposed activity-ratio distribution e �μ(σ ) [14]. Our interest in the influence of polydispersity on the hard-sphere fluid–solid transition led us to consider a quadratic distribution of chemical potential differences [9] β �μ(σ ) = −(σ − σ0)/2ν. (2) In the limit ν → 0, the pure monodisperse σ0-phase can be recovered. For small ν the mixture is ideal and the composition will be Gaussian with the peak near σ0. However, as ν increases and the distribution of diameters broadens, the mixture departs from ideal behaviour and the composition can no longer be expected to coincide with the imposed activity distribution. In this instance we find that the original form taken for the quadratic loses its relevance, and it becomes appropriate to write a more general form for the (stillquadratic) distribution β �μ(σ ) = c2σ 2 + c1σ + c0. (3) Obviously, equation (2) is recovered by identifying c2 = −1/2ν, c1 = σ0/ν and c0 = −σ 2 0 /2ν. The choice of equation (3) converts Y from a functional of �μ(σ ) to a function of σ0 and ν, and the fundamental thermodynamic equation can be written as d(βY ) = H dβ + βV dP + βμ(σ0) dN − (Nm1/ν) dσ0 − (Nm2/2ν) dν (4) Freezing in polydisperse colloidal suspensions 9629 where m1 and m2 are, respectively, the first and second moments of the composition about σ0. In a semi-grand Monte Carlo simulation, particles sample diameters in addition to the usual sampling of positions within the simulation box. Diameters are sampled by selecting a particle at random, changing its diameter by a small amount, and accepting with probability in accord with the Metropolis algorithm. Details may be found in [14]. Evaluation of the hard-sphere fluid–solid coexistence line as a function of polydispersity can be done by application of the Gibbs–Duhem integration method recently developed by one of the authors [12]. The governing differential equation in the present instance is in the space of P and ν [9] � dP dν � = �m2 2ν2β �v . (5) We integrate along the coexistence line in the (P, ν) plane by measuring the second moment m2 of the composition distribution and the molar volume v in both phases and applying the predictor–corrector scheme described elsewhere [12]. When integrating in ν from the monodisperse limit the method is helped by knowledge of the slope of the coexistence line at the initial point of the integration. This is not measured in a straightforward fashion, as the right-hand side of equation (5) is the ratio of quantities that vanish in the limit of a monodisperse system. Nevertheless, the ratio can be measured via tabulation of the distribution of particle nearest-neighbour distances in the monodisperse system [9]. Introduction of scaled coordinates and diameter allows (near-) analytical volume integration which greatly enhances the accuracy of the integration [9]. We performed a Gibbs–Duhem integration in the (P, ν) plane starting with an fcc solid and a fluid both consisting of 256 particles at the monodisperse hard-sphere freezing pressure. We evaluated the coexistence line by gradually increasing ν from one simulation to the next while integrating equation (5) to determine the pressure. A surprising result from the integration is a naturally occurring maximum in the polydispersity parameter ν, which arose as the P–ν curve bent back on itself, followed by a divergence of the pressure while ν continued back toward zero. This divergence of the pressure is perplexing until one realizes that the system is not returning to the monodisperse limit with the decrease in ν. Instead in this ‘infinite-pressure’ limit one finds that the particle diameters are vanishing as the integration proceeds, and that all quantities need to be rescaled by a length characteristic of the particle sizes in order for everything to remain well behaved. The parameter σ0 loses all significance as a characteristic length in this limit. In its place we find that the group ν/σ0 arises as the only important characteristic length (other than that associated with the pressure, (βP)1/3). All quantities, when scaled by this group, become finite in the ‘infinitepressure’ limit. In equation (3), this group is identified as the reciprocal of c1. Choosing it as our unit of length is equivalent to taking c1 ≡ 1. Importantly, when c2 is scaled in accordance with this, we have c2/c 1 = −ν/2σ 2 0 , which vanishes when ν → 0 in the infinite-pressure limit. Thus the terminal polydispersity that we identified in [9] corresponds to a system that has a chemical potential distribution described by equation (3), with c1 = 1 and c2 = 0. This system has a chemical potential distribution that rises unbounded with increasing particle diameter σ . This means that the imposed distribution is encouraging the particles to adopt larger diameters. Opposing this is the imposed pressure. The average volume and particle size distribution adopted by the system will depend only on the ratio βP/c3 1. Our calculations show that the coexistence pressure at the terminal polydispersity βP/c3 1 is approximately 0.000 15, while the average diameter σ/c1 adopted by the particles in each 9630 P G Bolhuis and D A Kofke Figure 1. The phase diagram in the plane of volume fraction and polydispersity. Coexisting phases are joined by tie-lines, which are not straight because the polydispersity is not an additive variable. The circles represent the end-points of the coexistence region for βPσ 3 0 → ∞, i.e., the terminal polydispersity. phase is of the order of 50. The distributions of particle diameters in the coexisting phases at the terminal polydispersity indicate that there is a significant degree of fractionation. The fluid adopts a wide distribution of smaller particles (on average), while the solid takes on a relatively narrow distribution of larger ones. It is convenient to choose a density variable η given in terms of the real volume fraction, because this is the quantity one measures in experiments. We define η = Nπ�σ 3�/(6V ); note that η = 0.7405 for monodisperse close-packed spheres. The phase diagram in the (η, s) plane is shown in figure 1. Because the polydispersity variable is not linearly additive, the tie-lines are curved. Coexisting fluid and solid hard spheres will mix to density– polydispersity values given by the lines (the converse is not necessarily true, as there are many compositions possible for a particular s, and an arbitrary one may not split exactly to the coexisting phases computed for the diagram). At the terminus the polydispersity is 5.7% for the solid and 11.8% for the fluid while the volume fractions are 0.588 and 0.547 for the solid and fluid respectively. We emphasize that the choice c1 = 1 made at the terminal polydispersity is done with no loss of generality in the results. As c1 and the pressure βP provide the only length scales, choosing c1 = 1 is like setting the diameters of a system of pure hard spheres equal to unity, σ = 1. Thus the only way to meaningfully perturb the system from this state is by introducing another non-zero term in the chemical potential distribution, equation (3). The least arbitrary choice is to work with non-zero c2. We will report on the influence of Freezing in polydisperse colloidal suspensions 9631 non-zero c2 in a future paper.